Interaction between domain walls in chiral p-wave superconductors 
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We calculate microscopically the interaction energy of domain walls separating degenerate ground 
states in a chiral p-wave superconductor. The interaction is mediated by the quasiparticles expe- 
riencing Andreev scattering at the domain walls. As a by-product, we derive a useful general 
expression for the free energy of an arbitrary nonuniform texture of the order parameter in terms 
of the quasiparticle scattering matrix. 
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I. INTRODUCTION 

Recent years have seen an increase in interest in topo- 
logical superconductors and superfluids from both exper- 
imental and theoretical facets,^ with one of the most 
studied examples being the chiral p-wave triplet state. 
The defining property of these systems is that, while 
the fermionic excitations in the bulk are fully gapped, 
there are gapless quasiparticles, which are protected by 
topology and are localized near inhomogeneities of the 
order parameter, such as sample boundaries, domain 
walls (DWs), and Abrikosov vortices. The chiral p-wave 
state in particular has received notable consideration be- 
cause of the gapless quasiparticle excitations (Majorana 
fermions) and non-Abelian winding statistics associated 
with half-quantum vortices, which are potentially useful 
as a route to quantum computing.^ 

The chiral p-wave triplet pairing, with order parame- 
ter proportional to kx ± iky, is experimentally realized 
in the superconducting state of Sr2Ru04 (Refs. U and 
0), as well as the A-phase of superfluid ^He (Ref. HI). 
The ground state is two-fold degenerate in the absence 
of an external magnetic field, and this gives rise to the 
possibility of superconducting or superfluid states with 
opposite chirality, separated by DWs, to form in differ- 
ent parts of the systemi^ii There is in fact experimen- 
tal evidence of the existence of DWs in Sr2Ru04 from 
Josephson measurements?^ and also in thin films of '^He- 
A from torsional oscillator measurements,^" In general, 
pairing states in other unconventional superconductors 
can exhibit discrete degeneracies of the ground state^ii 
which also leads to the possibility of DW formation in 
these systems. 

The formation of a DW costs gradient energy to the 
system due to the spatial variation of the order param- 
eter. Unlike ferromagnets, which break up into domains 
in order to minimize the net magnetic moment, in a neu- 
tral superfluid there is no analogous energetic rationale 
behind the formation of DWs. One possible explanation 
is that the DWs are spontaneously formed due to sample 
inhomogeneities during cooling across the phase transi- 
tion into the superfluid state. Alternatively, the creation 
of low-energy quasiparticles bound to the DW may com- 
pensate for the increase in gradient energy, which is par- 



ticularly effective in one-dimensional systems>ii 

The purpose of this paper is to develop a general micro- 
scopic formalism for calculating the interaction between 
superconducting DWs at arbitrary temperature. The 
structure of a single DW was investigated in Refs. 13, 
and|i2|. In general, the structure of superconducting DW 
textures can be studied using the Ginzburg-Landau (GL) 
formalism. It turns out that there are no stable two-DW 
solutions, and consequently, there must be some form 
of interaction between the DWs. Therefore, either an 
attraction between two DWs must cause an effective col- 
lapse of the DWs to a single domain; or a repulsion be- 
tween them pushes one of the DWs to infinity, leading 
to the effective formation of just two domains. It is this 
interaction which has stimulated our current work. 

The paper is organized as follows: In Sec. |lTl we in- 
troduce the two-DW configuration, for which we will ul- 
timately determine the interaction energy. In Sec IIII[ 
we compute the quasiparticle spectrum of this texture 
in the semiclassical (Andreev) approximation. Then, us- 
ing the transfer matrix method, we relate the interaction 
energy to the scattering matrix of the Bogoliubov quasi- 
particles in Sec IIVI Finally, we analytically calculate the 
interaction energy between DWs in the limit of large DW 
separation. Throughout this paper, we use the units in 
which h = ks = ^■ 



II. 



THE MODEL 



We consider a two-dimensional chiral p-wave neutral 
superfluid. Any external fields and disorder are ne- 
glected and an isotropic Fermi surface is assumed. The 
order parameter of a triplet fermionic superfluid or su- 
perconductor is a 2 X 2 spin matrix which has the form 
A(fe,r) = i(T2ard{k,r), where & are the Pauli matrices 
and d is the spin vector. For unitary states the latter de- 
fines the normal to the plane in which fermions paired at 
(fe, —k) are equal spin paired.^ In our case, d has only z- 
component, and its momentum dependence is given byi^ 



d = 



(1) 



where rji and 772 are components of a complex order pa- 
rameter vector ?7 and /cf is the Fermi wave vector. 
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We focus on planar superconducting textures describ- 
ing one or more DWs perpendicular to the x axis, there- 
fore only x-dependence is retained in tj. The DWs sepa- 
rate states of opposite chirality, hence the order parame- 
ter alternates between kx -I- iky and k^ — iky states. The 
spatial dependence of r] can be studied using, e.g. the 
GL formalism, see Appendix. There is no exact analyt- 
ical solution for the DW structure, even in the case of 
a single DW, and a variety of approximations have been 
proposed in the literature (Refs.[10, andfl^. Most im- 
portant qualitative features of the DW textures can be il- 
lustrated using the constant-amplitude model introduced 
by Volovik and Gor'kov in Ref. l6|. In this model the or- 
der parameter has the form ri{x) = Ao(l, e~*'''*^^'')e''^'^^\ 
where (p is the common phase, and 7 is the relative phase, 
of the order parameter components. 

The phases cf) and 7 are not independent. Conservation 
of current requires that the transverse current is constant, 
and since it is fixed by external sources, one may set it to 
zero. This results in a linear relationship between the gra- 
dients of </) and 7. Thus the DW texture can be described 
in terms of a single variable - a spatially-dependent rela- 
tive phase 7(a;). Variational minimization of the GL free 
energy functional with respect to 7 leads to a sine-Gordon 
equation, whose simplest nontrivial solution correspond- 
ing to a single DW has a kink-like form, connecting the 
asymptotics 7(±cx)) = ±7r/2 and varying within a re- 
gion of thickness S,d (which has the meaning of the DW 
thickness). In general, different models give different ex- 
pressions for ri{x), but the the condition 7(±oo) = ±7r/2 
always holds without reference to a specific profile of the 
order parameter near the wall. Furthermore, the com- 
mon phase difference between the two domains is fixed by 
the condition of vanishing supercurrent across the DW, 
see Appendix for details. Thus, one can write the order 
parameter asymptotics far from the single DW as follows: 



J7(x) = Ao(l,i), 



at X 



-00, and 



77(x) = Aoe^>^(l,-i), 



(2) 



(3) 



at a; — >■ -|-oo. Here x is a parameter depending on the mi- 
croscopic details of the system, satisfying the condition 
< X < TT. One can make analytical progress by consid- 
ering the sharp DW model, in which case — >■ and the 
order parameter changes abruptly at x = between its 
asymptotic values. 

We now consider two DWs at a fixed separation L, 
with the first DW positioned at 2: = 0, and the second at 
X — L. Using a similar setup as in the single DW case, the 
chirality alternates between the three domains, and we 
analogously impose the constraint of vanishing supercur- 
rent along the x axis, which leads to a non-zero common 
phase difference between the domains. The outer left re- 
gion {x < 0), and the region on the far right {x > L), 
correspond to the kx+iky state, while the middle domain 




FIG. 1: Alternating chirality states in the two-DW model. 
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FIG. 2: The relative phase between the order parameter com- 
ponents for two sharp DWs with fixed separation L. 



(0 < X < L) corresponds to the kx — iky state, as shown 
in Fig. m 

As in the single DW case, we focus on the sharp DW 
model to obtain an analytical solution for the interac- 
tion energy of the two DWs. Then the order parameter 
for both of the outer domains is given by the expression 
in Eq. and a non-zero global phase factor appears 
in the order parameter of the middle domain, which is 
given by the expression in Eq. ([3]). In accordance with 
the sharp DW model, 7(0;) changes abruptly between its 
asymptotic values in the three domains, as illustrated in 
Fig.H 



III. QUASIPARTICLE SPECTRUM 

Since we consider a neutral superfluid, interaction be- 
tween the DWs can only be due to their effect on the 
Bogoliubov fermionic quasiparticles. The quasiparticle 
spectrum for a nonuniform superconductor is determined 
by the Bogoliubov-de Gennes (BdG) equations, with the 
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4x4 BdG Hamiltonian given by 



n = 



At -I 



kp/m*, and vp^x = vpcosd. The DW order parameter 
at given kp has the form 

(4) Afc,(x)=ryi(a;)cos0 + ?72(x)sin0 = |Afc,(a;)|e**(-', (8) 



where ^ — ^{k)(jQ. Here ctq is the 2x2 identity matrix, 
and we assume that ^(fc) — {k^ — kp)/2m* , where m* is 
the effective mass. For the chiral p-wave states, we have 
A — dzCFi, where cx i iky. This form for the gap 
function allows the BdG Hamiltonian in Eq. (U), which 
operates on a four-component quasiparticle wavefunc- 
tion, to be written as a direct sum of two identical 2x2 
matrices, denoted by HbcIG- The four-component wave- 
function is thus decoupled into two two-component wave- 
functions which satisfy Hsda'^a = E"^^, where cr = ±, 
and HbcIG is given by 



H 



BdG 



4 



(5) 



We should point out that a does not denote the spin pro- 
jection of the two-component wavefunction; in fact, the 
components of both spinors have mixed spin projections. 
From this point on, we may drop the label cr, whose only 
effect is to double the degrees of freedom. 

For a DW parallel to the y-axis, where the order 
parameter depends only upon x, the y-dependence of 
the quasiparticle wavefunction is equivalent to that of 
a free particle (i.e. a plane wave). It can be written as 
e''=»yv[/(a;), where ^'(a;) satisfies the two-component BdG 
equations for a given ky: 



2m* 



Aix) 
2m* 



(6) 



with kx = — iV:E, fco = 



fc^, and A(x) = dz{x) = 



mix)ikx/kF) + ri2{x){ky/^F), see Eq. 

The DW order parameter A(a;) varies slowly on the 
scale oil/kp. Consequently, we can apply the semiclassi- 
cal (Andreev) approximatioi*^^ and seek solutions of the 
form "^{x) = e^''^^ip{x), where ip = (u, u)'^ is a slowly 
varying "envelope" function with electron-like (u) and 
hole-like (w) components. Due to the circular symmetry 
of the Fermi surface in the xy plane, we have kx = ±fco 
for a given ky. Substituting ^(o;) of this form into Eq. ([6]) 
and neglecting terms containing higher-order gradients of 
ip we see that the envelope function satisfies the Andreev 
equations: 



-ivp^x^x Afep(x)' 
AL(a;) ivp^x^x 



ip = Eip. 



(7) 



The direction of semiclassical propagation of quasiparti- 
cles is defined by the Fermi wavevector kp = {kx,ky) = 
kp{co86,sm9), the Fermi velocity is given by vp = 



where |Afcy| is the magnitude of the semiclassical order 
parameter and $ is its phase. 

The asymptotics of A^^ (x) for the chiral p-wave state 
are fixed by Eqs. © and ([3]); however, different mod- 
els for the DW structure, see Sec. [TTl lead to different 
forms for the order parameter in the vicinity of the DW. 
In the single sharp DW model, the order parameter is 
uniform within the domains of opposing chirality, chang- 
ing abruptly at the boundary a; = 0. Thus we have 
Afc^(x) = A^e{-x) + A+9{x), where A_ = Aqc*^, 



A, 



Aoe^ 



and 9{x) is the Heaviside step func- 



tion. For two DWs, we have A^e{-x) + A+9{x)9{L - 
x)+ A_9{x~ L). 

There are two types of solutions supported by the 
Andreev equations ([7]): discrete bound states (Andreev 
bound states, or ABS's) for which \E\ < Aq, as well as 
a continuum of scattering states where \E\ > Aq. It will 
be shown below that all quantities of interest, including 
the interaction between the DWs and also the ABS spec- 
trum, can be expressed in terms of the properties of the 
scattering states, encoded in the scattering matrix S. 



A. Scattering matrix for two DWs 

The scattering matrix S relates the amplitudes of the 
incident wavefunctions to the outgoing amplitudes of the 
waves that are reflected/transmitted by the DW config- 
uration. To facilitate the calculation of S, we perform 
a gauge transformation on the wavefunctions ip to re- 
move the phase in the order parameter A^^ , see Eq. ([8]) , 
that appears in the off-diagonal elements of the Andreev 
Hamiltonian. For the two-DW setup described in SecHIl 
we can adopt a more convenient notation to describe the 
order parameter in each of the three domains: 



Afe^(x) 



Aoe 
A, 



oe" 



a; < 0, X > L 
< X < L 



(9) 



which can be further simplified to Akp{x) = Aoe**(^\ 
with $ = (^2 = X ~ ^ in the middle domain, and $ = 
ipi — 9 in the outer two domains. 

We denote the gauge-transformed wavefunctions by 
ip^x) and define ^p = Uip, where U is given hy U — 



The gauge-transformed wavefunctions sat- 



„i*(a:)CT3/2 

isfy 'ip{+oo) — ip{—oo), which follows from the condi- 
tion ip{-\-oo) = '0(— oo). The latter is consistent with 
the gap equation and the order parameter asymptotics 
Afcp(+oo) = Afc^(-oo). 

Continuity of the original wavefunctions at the bound- 
aries X = 0,L implies V'(+0) = 0) and ipiL -|- 0) = 
ijj{L — 0); however, removing the phase from the order 
parameter causes the gauge-transformed wavefunctions 
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tjj to suffer a phase discontinuity at the domain bound- 
aries. In fact, one can easily verify that ip satisfy the 
foUowing conditions: 



where 



^(+0) 



-tS 



^(-0), 



(10) 



where S — [ip2 — </'i)/2 = x/2 — 6. 

Ahhough we remove the phase $ from the off-diagonal 
terms in the Andreev Hamiltonian, its derivative $' ap- 
pears in the diagonal elements after the gauge transfor- 
mation, so that -0 satisfies the equation 



-ivF,xCrzVx + ]^vf,x^' {x)aQ + Ao(7i 



V' = i?V, (11) 



with the same energy eigenvalues as the original wave- 
functions. Due to the delta-function singularities of $' 
at a:; = 0, L, we apply the gauge transformation sepa- 
rately in each region (i.e. at all a; 7^ 0, L, where $' = 0). 
In this way we obtain an equation in each domain which 
has the form of the Andreev equation in a uniform super- 
conductor, and as such can easily be solved. The solution 
must satisfy the "twisted" matching conditions given by 

Eq. (nni). 

We focus on the scattering states and from now on 
drop the tilde on the gauge-transformed wavefunctions. 
For the continuum of scattering states, the quasiparticle 
wavefunctions are linear combinations of plane waves: 



'lj}{x) 



(12) 



where u± = Ao/(£' T gwF,x), with q = ^/W^^^/\vf,x\ 
> 0. The subscript on the amplitudes in Eq. ((T^ corre- 
sponds to the direction of quasiparticle propagation (i.e. 
left or right). We also introduce a superscript on the am- 
plitudes of the wavefunctions in the outer two domains 
(x < 0, x > L) to identify each particular region. Let 
the "— " superscript denote the region a; < 0, and the 
"-h" superscript correspond to the x > L region. Then 
the amplitudes of the waves incident on the DW config- 
uration are given by ■* and a''^^ , while the outgoing 



(reflected and transmitted) waves have amplitudes A\_ 
and ylfT''. We define the scattering matrix S as follows: 



aV\ 



(13) 



The scattering matrix is calculated by using the match- 
ing conditions given by Eq. (flUl) to eliminate the wave 
amplitudes in the region < a; < i and relate the ampli- 
tudes of incident waves to those of the outgoing waves. 
The final result has the following form: 

5*11 ~ S22 = -p, 

R- R+ 
^12 - -p-, ^12 - -p-, 



P = 1 



^0 



2iqL 



l)sin2 6, 



R± = ig±l) {icosS±gsmS) [e^^'i^ - l)sm6, 

and g = E/qvp.x- 

To conclude this subsection we note that, by consider- 
ing scattering from the left and right separately, one can 
relate the scattering matrix to the reflection and trans- 
mission coefficients of the Bogoliubov quasiparticles in 
the presence of the order parameter texture. If, for exam- 
ple, there is a wave incident on the DW located at a; = 
from the left, then we can set the incident amplitude 
A^"-' = 1, and we must have A*!^' = 0, A^^^ = and 

^ — tl, where and represent the left-incident 
transmission and reflection coefficients, respectively. The 
right-incident transmission and reflection coefficients Ir 
and rn can be introduced in a similar way, by consider- 
ing right-incident scattering on the DW at a: = L. Then 
one can relate the scattering matrix entries to the trans- 
mission and reflection coefficients as follows: = Sn, 
tl = 5*21, tR = S22, and vr = Si2- We should point 
out that since we have applied the gauge transformation 
to the quasiparticle wavefunction before calculating the 
scattering matrix, the reflection and transmission coef- 
flcients obtained by this method are not equivalent to 
the ones that would be obtained from the direct Andreev 
calculation (prior to the gauge transformation). 



B. Bound state spectrum 

The ABS energy for the single sharp DW model has 
been previously calculated,— and it has the following 
form: 



Eo{e) = AQs{e) cos (^6 



(15) 



with s{6) — sgn [sin (6* — x/2) cos 6*]. This expression is 
valid for an arbitrary phase difference across the DW and 
it should be noted that, in general, the ABS energy is not 
a continuous function oiO. There are certain directions of 
semiclassical propagation at which discontinuities occur: 
at 6 = ±7r/2, corresponding to a "grazing trajectory" 
where the quasiparticles move parallel to the DW (in this 
case the Andreev approximation is actually not applica- 
ble); and also at = x/2 and 9 — x/2 + tt, in which case 
the quasiparticles do not "see" the DW, since A_|. = A_ . 
The ABS energies for x = 0, tt were calculated in Ref . [TgI. 

We proceed now with the calculation of the bound 
state energies for two DWs. The ABS energies are 
obtained from the poles in the scattering matrix en- 
tries, see Eq. after analytical continuation to the 
real energy axis within the interval \E\ < Aq. Defln- 
ing the dimensionless energy e — E/Aq, we rewrite 
(14) q = (Ao/wF)Ve^ — 1/| cos6\. For the bound state ener- 
gies with |e| < 1, we have to choose the correct branch of 
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q before proceeding further. To this end, we consider the 
function w{z) — \J — \^ which has two branch points: 
one at z = 1, and the other at z — —1. We choose the 
branch to ensure that w(z) is real and w(z) > if z is real 
and \z\ > 1, in accordance with the definition of q. One 
can select the branch cuts to run parallel to the imagi- 
nary axis, from ±Ao to ±Ao =F and then the correct 
choice \bw{z) = i\J\ — for z along the real axis within 
the interval [—1,1]- 

To simplify the denominator P in Eq. ([T4| . we define 



— ,„>^ — 



5(e) 



'I, where L = Lj^ is the 



dimensionless distance between the DWs and ^ = vp/ 
is the correlation length. Then from the poles in the 
scattering matrix entries we obtain the following equation 
for the bound state energies: 



— cos^ 5 — (5(e) sin^ 5 = 0. 



(16) 



To make analytical progress, we can consider this last 
equation for small values of a, which physically corre- 
sponds to large DW separation. Expanding the equation 
in powers of a up to linear order, we obtain the result 
e = ±1 cos5|[l -I- (atan^ 5)/2]. The two DWs become de- 
coupled at a — > 0, which occurs for either the grazing 
trajectory 9 = ±7r/2, or for L 1. In this case, we 
recover the same shape for the ABS energy as that as- 
sociated with the single DW configuration, see Eq. ([T5|). 
but the dependence on s{9) has been lost. The bound 
state energies for the DW located at a; = correspond to 
Eq{9), while those for the anti-wall at x — L correspond 
to —Eq{9), and consequently, our energy curves for the 
two DWs have two branches, as shown in the top panel 
of Fig. 13] in the case of x = 

We solve Eq. numerically to obtain a profile for the 
bound state energies in the general case. To illustrate the 
effect of L on the spectrum, we present the result for x = 
TT in Fig. |3l As the distance between the DWs decreases, 
the bound states localized near them become hybridized, 
which leads to the splitting of the energy branches. 



IV. INTERACTION BETWEEN DOMAIN 
WALLS 

We are now in a position to evaluate the interaction 
energy between the DWs. We begin with the expression 
for the free energy of an arbitrary nonuniform supercon- 
ducting texture in terms of the Fredholm determinant of 
the BdG Hamiltonian. For a system with two-component 
order parameter in zero magnetic field, it has the form^i 



J" = 



-T^lnDet'*'^""^^'^^ 



1 

'v 



{\rii\' + \V2\')d\ 



(17) 



where is the total free energy of the system measured 
with respect to the normal state with A(x) = 0, HbcIG is 



E/A„ 



FIG. 3: Bound state energy for the two sharp DW model with 
X = TT, for varying dimensionless DW separation L — L/^. 
From the top: i = 5, 1, 0.5. 



the BdG Hamiltonian defined by Eqs. (IS|) and is 
the normal-state BdG Hamiltonian, ojn = [In -f- 1)ttT is 
the fermionic Matsubara frequency, and V is the coupling 
constant in the chiral p-wave channel. 

We introduce the free energy difference, ^ between 
nonuniform and uniform superconducting states, where 
the nonuniform state has two DWs. From Eq. (fT7|) . it 
follows that 



5T = -T^lnDet 



ILOn — HBdG 



IUJ„ 



H 



(0) 



1 



BdG 

J {\V\'-W'^\')d\ (18) 



where 77*^°^ and H^^q denote the order parameter and 
the BdG Hamiltonian corresponding to the uniform chi- 
ral state. The expression ^TE\\ depends on the separation 
between the DWs. Since the DWs are decoupled at infi- 
nite separation, SJ^{L — >• 00) gives the self-energy of two 
DWs, measured with respect to the uniform supercon- 
ducting state. The interaction energy between the two 
DWs Tint is simply the difference between the free energy 
at arbitrary DW separation L and the self-energy of the 
two-DW configuration, i.e. J-mt = SJ-{L) — SJ-{L — 00). 

For the sharp two-DW model introduced in Sec. [TTl we 
have |?7p = 2Aq in each of the three domains. Also, 
1^(0) |2 _ 2Aq, and consequently the second term on the 
right-hand side of Eq. ([T5)) vanishes, and the interaction 
energy takes the form 



where 



J"(L) = -T^lnDet 



H 



BdG 



ILUri 



H 



(0) 
BdG 



(19) 



(20) 
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The logarithm of each of the Fredholm determinants in 
Eq. can be written as follows: 



In Dct 



- H 



BdG 



H 



(0) 

BdG J 



e: 



(0) 



(21) 



where « is a set of quantum numbers labelling the eigen- 
states of the 2x2 BdG Hamiltonian, see Eq. ([6]) , at given 
DW separation L, Ei are the corresponding eigenvalues, 
and Ef'' are the eigenvalues for the uniform chiral state. 

The sum over the BdG spectrum in Eq. (|21l) can be ex- 
pressed in the semiclassical approximation as a sum over 
the eigenvalues of the Andreev Hamiltonian Ha, defined 
by Eq. ([7]), as follows: 



2T:NFiy 



\VF.x 



(22) 



where Np — m/27r is the density of states at the Fermi 
level per one spin projection in two dimensions, £y is the 
length of the DW, kp defines the direction of semiclas- 
sical propagation of the quasiparticles, and j labels the 
eigenstates of the Andreev Hamiltonian at given kp. Re- 
call from section Sec. IHI Al that the quasiparticle wave- 
functions satisfy ^(-|-oo) = ■0(— oo), so we have appro- 
priately placed our system in a box of dimensions £x = i 
and iy and imposed the periodic boundary conditions. 

It follows from Eqs. ([Ill [201 1211 1221) that the interaction 
energy per unit DW length is given by: 



^Vp,.j:\\n- 



27r' 



L'(iw„; oo) ' 



(23) 

with D{z) = Y{^{z ~ Ej)/{z - Ef''), where Ej are the 
eigenvalues of the Andreev Hamiltonian at given kp 
for given DW separation L, and are the eigen- 

values of the Andreev Hamiltonian in the uniform chi- 
ral state. Note that the gauge transformation intro- 
duced in Sec. IHI Al leaves the eigenvalues of Ha unaf- 
fected. After the transformation, the Andreev Hamil- 
tonian can be written as Ha — H^a' ^ where 
H^2^ = —ivp^x<^^'^ X + is the Andreev Hamiltonian 

for the uniform chiral state which now has A^^ (x) = Aq, 
and 5H = vp^^^' {x)crQ /2 is a localized perturbation, cf. 
Eq. pT|) . Adiabatically switching on the perturbation, 
there is a one-to-one correspondence between the eigen- 
values of Ha and H^a ■ 

Next we introduce the 2x2 transfer matrix M{x\ E) 
which acts as an x-evolution operator for the quasi- 
particle wavefunctions at given energy E: ip{x) = 
M{x; E)ip{—£/2). The transfer matrix satisfies the fol- 
lowing conditions: 

iHA-E)M{x;E)^0, M {~e/2; E) ^ ao, (24) 

which hold for arbitrary values of E. It follows from 
the periodic boundary conditions that the quasiparti- 
cle wavefunctions satisfy [ao — M{£/2; E)]^{—£/2) = 



0. This quantization condition leads to the charac- 
teristic equation for the eigenvalues of Ha, given by 



det 



ao - M {£/2; E) 



0, where det (• • • ) is a 2 x 2 

determinant. We also define another transfer matrix 
M(°)(a;;-E), which satisfies the same conditions (|24[) . but 
for the uniform-state Hamiltonian h'^^^ . 

From here we can introduce a new quantity d{z), which 
is defined by the following expression: 



d{z) 



det ao-M{£/2;z) 



det 



ao- Ah{£/2;z) 



(25) 



Both d{z) and the Fredholm determinant D{z) have ze- 
ros at z = Ej, as well as poles sX, z = E, 



(0) 



For 



\z\ — > CXI, which physically corresponds to large values 
of E, the quasiparticles are not affected by the super- 
conducting order parameter. Consequently, Ha h'^a' , 
and D{z),d{z) — > 1. Due to these properties, we obtain: 



D{z) = d{z), 



(26) 



see, e.g. Ref. [l8| for review. 

In the following subsection we use the transfer ma- 
trix method, in particular Eq. (I^B)) . to relate the DW 
interaction energy to the scattering matrix entries. Sub- 
sequently, we evaluate the sum over the Matsubara fre- 
quencies and the integral over semiclassical directions of 
propagation in Eq. to obtain an analytical expres- 
sion for the interaction energy in the limit of large DW 
separation. 



A. Calculation of the Fredholm determinant 

To facilitate the calculation of the Fredholm determi- 
nant at imaginary (Matsubara) energies, we first define 
a matrix f , which relates the amplitudes of the gauge- 
transformed quasiparticle wavefunctions on the left-hand 
side of the DW configuration to those on the right-hand 
side, see Sec. IHI Al as follows: 



(+) 



i-y 
+ 

(-) 



(27) 



It can, therefore, be expressed in the following way: 
1 fdetS Si2\ 
^ " S22 [-S21 1 / ' 



(28) 



where S is the scattering matrix defined by Eq. (jl3p . 

We introduce a shorthand notation for the transfer ma- 
trix from -£/2 to +£/2: M{£/2; z) = m. Using Eq. ^ 
and the wavefunctions defined by Eq. p2)) . we find that 



V+tV_ , where 



V± 



±iqe/2 
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and u± and q are defined in Sec. IIII Al We introduce 
a similar notation for the uniform-state transfer matrix 
M(°)(£/2; z) — mo, and since in the absence of DWs f = 
S = aQ, it immediately follows that mo = V^VZ^- 

We can now rewrite the expression for the Fredholm 
determinant, see Eqs. (|25]) and (l26l) . as 



Diz) 



det ( (To - V+fVl'^ 



det CTo - 



(29) 



and after multiplying the matrices we find that the nu- 
merator in the last equation takes the form 



det ((To - m) = 1 + det f - (rne'^'' + T22e 



(30) 



In the uniform superconducting state, this reduces to 
det (iTo — mo) = 2(1 — cosqi). 

To calculate the interaction energy, we must evaluate 
D{z) at discrete imaginary points z = iujn, see Eq. p3|. 
and since q in Eq. ([50]) is only defined for real values of E 
such that \E\ > Aq, we must analytically continue q in 
the complex energy plane to the imaginary energy axis. 
Using the same procedure as in Sec. IIII Bl we find that 
the appropriate expression for q is q{E — iuJn) = ix, 
with K — -y/a;2 -I- 1^^/\vf,x\- Therefore the exponential 
terms in Eq. ([50)) take the form e**^^ = e^'^^, and in the 
thermodynamic limit ^ — oo, we keep only the last term 
in this equation as the others are small in comparison. In 
this limit, the Fredholm determinant in Eq. (|29p becomes 



The overall sign of the interaction energy is negative, so 
we see that the DW interaction mediated by the Andreev 
scattering of quasiparticles is attractive at all tempera- 
tures, which means there is an effective collapse of the 
walls to a single uniform domain. Qualitatively it is ev- 
ident from Eq. p3|) that the attraction is exponentially 
weak in the limit of large separation. To make analytical 
progress, we focus on the case of zero temperature and 
X = 0. The results for other values of the phase difference 
are expected to be qualitatively similar. 

At zero temperature, the summation over the discrete 
Matsubara frequencies aj„ becomes an integral over a con- 
tinuous variable uj: T'}2,^[- •■)—>■ /(• • • )da;/27r. In the 
limit of large DW separation, one can expand the loga- 
rithm in Eq. (|33p, then the interaction energy takes the 
form: 



27r 



2tt 



I COS 6*1 



+ 00 



sm 



-t- COS^ I 



■ exp 



2L^/^ 



I COS 6*1 



.(34) 



where Co = w/Aq, and the dimensionless distance L was 
introduced in Sec. IIII Bl In the limit L ^ 1, one can 
further neglect the w-dependence of the pre-exponential 
factor in the integral over Co, and evaluate this integral 
by the steepest descent method. In this way, we can 
represent Eq. (|34[) in the following form: 

2Nf^ovf 



(35) 



D{z = iuJn) = T22{iuJn) 



522 (iw„) 



(31) 



where we have used the relation between f and S given 
in Eq. ([28|) . Thus the calculation of the Fredholm deter- 
minant has been reduced to finding the properties of the 
scattering states. 

The expression ([51]) is applicable to any planar super- 
conducting texture. In the case of two sharp DWs, using 
Eq. ([T4)) . we obtain: 



1 



Al sin^ S 
i^l + Al 



1 



-2>{L\ 



(32) 



where 5 was defined in Sec. IIII Al 



B. Interaction energy 



From Eqs. and (|32|), it follows that the in- 

teraction energy per unit DW length at given DW sepa- 
ration has the following form: 



F,nt = -vfNfTY] / 

„ "'0 



d6'|cos6'| 



X In 



A^ sin-^ 5 



Aq cos^ 5 



(33) 



where 



I{L) 



7r/2 



sm 



,-2L/|c 
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VI cos 6*1 



(36) 



and we have invoked the symmetry of the angular integral 
to reduce the region of integration. 

We make a change of variable p = 1/ cos 9 in Eq. 
and obtain: / = dp{p^ - l)i/2p-5/2g-2Lp_ ^j^jg j^^g^ 
integral can be expressed in terms of the modified Bessel 
functions of the second kind, and in the limit i 3> 1, 

it has the form / = (10/3)\/ 7rLe~^^. Using this result 
in Eq. psp . and restoring the dimensional quantities, we 
arrive at our final expression for the interaction energy 
per unit DW length: 

^ (37) 



Fi, 



20 A 

-—NfAqvf exp 



VF 



Thus, the interaction between the DWs is attractive and, 
as expected, it is exponentially weak in the limit of large 
separation between the walls. 



V. CONCLUSIONS 

We studied the interaction between two DWs separat- 
ing states of opposite chirality in a p-wave superconduc- 
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tor, and found that it is attractive for arbitrary DW sep- 
aration, at all temperatures. Furthermore, we found that 
the interaction energy is exponentially weak for large sep- 
aration between the DWs. We used the transfer matrix 
method to relate the interaction energy of the DWs to the 
scattering matrix of the Bogoliubov quasiparticles, and 
the latter was calculated in the semiclassical (Andreev) 
approximation. 

The transfer matrix approach developed in this paper 
has a more general validity and can be applied to any su- 
perconducting texture. The free energy can be expressed 
in the same form as Eq. ()18|) in terms of the Fredholm 
determinant of the BdG, or Andreev, Hamiltonian; how- 
ever, the quasiparticle scattering matrix will be sensitive 
to the order parameter configuration. It would be inter- 
esting to use this method to characterize the interaction 
between DWs in a variety of other unconventional su- 
perconductors with discrete degeneracies of the ground 
states. 

Our results are immediately applicable to the neutral 
case. In a charged superconductor, there will be an- 
other contribution to the DW interaction coming from 
the Meissner currents and associated magnetic fields. In- 
vestigation of these effects on the interaction energy is 
left for the future work. 
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Appendix A: GL Description of a Domain Wall 

To gain insight into the structure of the DW, we use 
the GL free energy functional. In a nonuniform neutral 
superfluid it is a sum of uniform {Fu) and gradient (Fg) 
energy densities. For the superconducting state with two- 
component order parameter r) — (?/i,?72), and isotropic 
Fermi surface, we have 



and 



F^ = a\r]\^ + P,\ri\^ + P2\VV\ 



(Al) 



(A2) 



where the Einstein summation convention is assumed. 
We find that the minimum in the uniform free energy 
density corresponds to the chiral states t] = Ao(l,±i) 
with Aq = ■\/|q;|/4/?i, for the case Pi, > 0. 

The order parameter of a planar DW configuration can 
be expressed in the following form: 



771(2;) -Ao/i(x)e^^(-), 
772(2;) -Ao/2(x)e*<^(-)-^^(-), 



(A3) 



where /i^2 are the dimensionless amplitudes of the order 
parameter components, whose asymptotics far from the 
DWs are given by /i^2 = 1- 

The nonzero phase difference x between domains of 
opposing chirality emerges from the condition of van- 
ishing supercurrent across the DW. It follows from 
Eq. (|A2[) that the supercurrent is given by the expres- 
sion = 2 Im {Ki-q*y .i-jj + K2T1*V j-qj + ^377* V^Ty^). Us- 
ing Eq. (jASp , we find that the transverse current has the 
form J, = 2A2(Ki23/i' +ifi/|)(V.0) - 2ifiA2/2(V,7), 
with K123. = Ki+ K2 + K3. 

In our work we focus on the static case, so conservation 
of current requires jx — const. One can set jx = 0, which 
yields a linear relation between the gradients of (p and 7. 
This relation can be used to eliminate the common phase 
from the gradient energy and obtain: 



F, 



F„ = 



aAUf! + fi)+PiAtif!+fir 

-K/32A4(/4 + /4+ 2/2/2 cos 27), 

i^l23A2(V,/l)2 + i^iA2(V,/2)^ 
KiKi23f?fi 



(A4) 



Ki23ft + Kif, 



Variational minimization of these expressions with re- 
spect to /1.2 and 7 gives rise to three coupled nonlin- 
ear differential equations. Using the solutions to these 
equations, one can compute the common phase difference 
between arbitrary points xi and X2'. 



(^2 



(^1 



Ki23ft + Kifi 



{Vxl)dx. (A5) 



We see that whenever there is a gradient of the rela- 
tive phase 7, the value of the common phase difference is 
nonzero, and it is evidently sensitive to the microscopic 
details of the system. While we have imposed the con- 
dition of zero transverse current to obtain Eq. (jASp . the 
current along the DW is nonzero. 

Obtaining an exact analytical solution for the DW 
structure is not possible due to the complexity of the dif- 
ferential equations obtained from variational minimiza- 
tion of the expressions in Eq. (jA4p . However, one can 
make analytical progress by considering the constant- 
amplitude modeljS in which case fi,2{x) = 1 for all 
X. Then it follows from Eq. (|A4p that the total free 
energy density is given by the expression F — F^ + 
Fg = (•••) + KAliVxl)^ + 2/32A4cos27, where K ^ 
K1K123/ (K123 + Ki). The first term in F contains 7- 
independent contributions, and it is immediately clear 
that the equation for 7 has a sine-Gordon form: 



ifV^7 + 2/?2A2 sin 27 = 0. 



(A6) 



The simplest nontrivial solution to this equation is 
a kink-like solution given by sin 7(2;) = tanh{x/^d)- 
This corresponds to a single DW. The parameter = 



K/4:(32Aq has the physical meaning of the DW thick- 
ness and is of the order of the GL correlation length. 
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Using this expression for 7(x) in Eq. (IA5[) we obtain: 
X = 0(+oo) — 4>{—oo) = ttKi/ {Ki23 + Ki). In the weak 
couphng model, Ki ~ K2 — K3 (Ref. [13), and x — t^/^- 
There are no stable two-kink (or two-DW) solutions 
to Eq. (jA6l) : the other nontrivial solution corresponds 
to a periodic lattice of DWs. This can be understood 
by considering a simple pendulum with 27 — >■ (the 
angular displacement of the pendulum), and x — )■ i (a 
time coordinate). Consider first the one- kink solution, 
where initially t — >■ —00, and we have = — tt. After 
a sufficient amount of time has elapsed, the pendulum 
has just enough energy to complete one full revolution, 
approaching an angular displacement of = -I-tt. There 
can be no two- kink solutions because if the pendulum has 



enough energy to surpass the limit & = +Tr and complete 
one more revolution, it will have enough energy to do 
this an infinite number of times, which corresponds to a 
periodic arrangement of DWs. 

The fact that there are no stable two-DW solutions 
implies there must be some form of interaction between 
the walls. If this interaction is attractive, it will cause 
a collapse of the DWs to a single domain. If it is re- 
pulsive, then one of the DWs will be pushed to infinity, 
leaving just a single DW separating two domains of op- 
posite chirality. The periodic solution for "f{x) can also 
be understood in terms of this interaction, as a mutual 
attraction or repulsion between neighbouring DWs could 
potentially lead to a stable periodic configuration. 
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